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ABSTRACT 

We describe the implementation of a module for the Athena magnetohydrodynamics (MHD) code 
which solves the time- independent, multi-frequency radiative transfer (RT) equation on multidimen- 
sional Cartesian simulation domains, including scattering and non-LTE effects. The module is based 
on well- known and well-tested algorithms developed for modeling stellar atmospheres, including the 
method of short characteristics to solve the RT equation, accelerated Lambda iteration to handle scat- 
tering and non-LTE effects, and parallelization via domain decomposition. The module serves several 
purposes: it can be used to generate spectra and images, to compute a variable Eddington tensor 
(VET) for full radiation MHD simulations, and to calculate the heating and cooling source terms 
in the MHD equations in flows where radiation pressure is small compared with gas pressure. For 
the latter case, the module is combined with the standard MHD integrators using operator-splitting: 
we describe this approach in detail, including a new constraint on the time step for stability due 
to radiation diffusion modes. Implementation of the VET method for radiation pressure dominated 
flows is described in a companion paper. We present results from a suite of test problems for both 
the RT solver itself, and for dynamical problems that include radiative heating and cooling. These 
tests demonstrate that the radiative transfer solution is accurate, and confirm that the operator split 
method is stable, convergent, and efficient for problems of interest. We demonstrate there is no need 
to adopt ad-hoc assumptions of questionable accuracy to solve RT problems in concert with MHD: 
the computational cost for our general-purpose module for simple (e.g. LTE grey) problems can be 
comparable to or less than a single timestep of Athena's MHD integrators, and only few times more 
expensive than that for more general (non-LTE) problems. 

Subject headings: (magnetohydrodynamics:) MHD methods: numerical radiative transfer 
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1. INTRODUCTION 

Radiation is of fundamental importance for the ther- 
modynamics of most astrophysical systems. It can be 
the dominant source of heating and cooling of astrophys- 
ical plasmas. Even in those systems where it plays a 
minor role in energy transport, it is the dominant mech- 
anism through which we perceive and explore the uni- 
verse. Nevertheless, it has often proven difficult to di- 
rectly model the effects of radiation accurately in modern 
multidimensional astrophysical (magneto)hydrodynamic 
(MHD) codes due to both computational expense and 
conceptual complexity. 

Most approaches to adding radiative transfer to dy- 
namical simulations are based on adopting restrictive 
assumptions or approximations. For example, often 
the flow is assumed to be optically thin to radiation 
everywhere and for all time, or the radiation field 
is assumed to originate in a small number of point 
sources, with the diffuse emission from scattered or rera- 
diated photons ig nored (such as in Cosmo l ogical reioniza- 
tion p roblems e.glAbel fc Wandeltl l2002t iMellema et all 



2006t iRiikhorst et all 120061: I WhaLen fc Normanl 120061: 



Reynolds et al.ll2009l: iFinlator et al.l l2009) 



For problems in which the diffuse emission cannot 
be ignored, the dynamics of the radiation field is often 
treated by solving the radiation moment equations using 
ad hoc closure prescriptions to handle the transition 
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from optically thick to opti cally thin regimes, such 
as flu x-limited diffusion (e.g. iLevermore fe Pomraninel 
119811 hereafter FLD). This includes applications such 
as accretion flows, star formation, neutrino transport in 
supernovae, stellar atmospheres and winds, cosmological 
reionization, and many others. Indeed, there is a large 
and growing list of astrophysical MHD codes that 
utilize FLD or a similar prescribed closure relation 



(including e 


.g.lTurner & Stonell2001l:lBruenn et al.ll2006t 


Haves et al. 


2006; Gonzalez et al.ll2007t iKrumholz et all 


2007al iGittines et all 


20081: ISwestv & Mvral 


20091: 
2TTTT1: 


Commcrcon et all 12011 


: Ivan der Hoist et al. 


Zhang et al. 2011). 



Numerical methods for directly solving rad iative trans- 
fer ( RT) have been implemented (e.g. IStone et all 



19921: Ivan Noort et all 120021: lHaves fc Normanl 120031 : 



Hubenv fc Burrowsll2007l ). but their application to astro- 
physical problems has been somewhat limited, especially 
in full 3D. A notable exception is the progress made in 
simulating the atmospheres of the Sun and other cool 
stars. In the solar physics community, multidimensional 
MHD simulations of convection with realistic RT have 
been performe d for decades w i th increasing sophistica- 
tion. fsee e.g. iNordlundl 119821: IStein fc Nordhmdlll99 



ILLS) 



Vogler et al iHeinemann et all 2007; H avek et al 



Encouraged by recent work modeling the departure 
of the radiation field from local thermodynamic equi- 
librium (LTE) due to the presence of electron scat- 
tering in th r ee-dim ensional MHD simulations (see e.g. 
I Havek et all l2010f ). we have implemented a general- 
purpose RT solver in Athena (jStone et alll2008l ). based 
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on the methods widely used in the stellar atmospheres 
community. Athena is a general purpose astrophysical 
MHD code, which is being actively developed and al- 
ready includes several modules for handling a variety of 
physical processes. Effectively, we have combined Athena 
with a modern stellar atmospheres code. In fact, Athena 
already has a RT module that computes the effects of 
ionization radiat ion from a single point source on the 
surrounding gas (jKrumholz et al . 2007b). However, this 
module is not well-suited for modeling the radiation from 
diffuse continuum emission. 

The addition of a RT solver to Athena enables three 
goals: (1) it can be used as a diagnostic tool to compute 
self-consistently spectra and images from time-dependent 
MHD simulations for direct comparison to astronomical 
observations; (2) it allows us to compute a variable Ed- 
dington tensor (VET) for the integratio n of the coupled 
MHD and radiation moment equations ([Sekora &: Stone! 
I201fi Jiang et al., submitted to ApJS, hereafter JSD12) 
for full radiation MHD simulations in regimes where both 
energy and momentum transport by photons is impor- 
tant; and (3) it allows us to compute the radiation source 
terms in the energy equations and directly couple them 
to the MHD integrator to compute the dynamics of flows 
where radiation pressure can be ignored. 

This paper focuses on describing our implementation 
of methods to solve the RT equation, and the coupling of 
the solver with the MHD integrator to compute the ra- 
diation source term in the energy equation. The compu- 
tation of the VET and solution of the radiation moment 
equations is described in JSD12. The plan of this work 
is as follows: In Section [2] we summarize the equations 
that are solved. In section [3] we describe the detailed 
implementation of our solver and the iterative methods 
used model deviations from LTE and handle certain (e.g. 
periodic) boundary conditions. In section |4] we describe 
how we compute the radiation source terms in the en- 
ergy equation and incorporate them into the MHD in- 
tegration. In Section [5] we present the results of several 
test problems not only to assess the accuracy of the RT 
solver, but also to evaluate the performance of the MHD 
integrator when the energy source terms are included. 
We summarize our results in Section [6l 

2. MHD EQUATIONS WITH RT 

In this work we solve the usual equations of compress- 
ible MHD, including the source term in the energy equa- 
tion to account for heating and cooling due to radiation. 
These source terms are computed directly from a formal 
solution of the time-independent RT equation. Thus, the 
basic equations are continuity 



f + V. ( pv)=0, 



momentum conservation 
d(pv) 



dt 



V • (pvv + T) = 0, 



the induction equation 
dB 
~dt 

and energy conservation 
dE 
~dt 



V x (v x B) = 0, 



V-(£v + T-v) =Q rad . 



(1) 



(2) 



(3) 



(4) 



In the above, p is the gas density, p, v is the fluid 
velocity, and B is the magnetic field. The total stress 
tensor T is defined as 



T = (p + S 2 /2)l-B T B, 
and E is the total (fluid) energy 



E 



P 



1 



-pv 



2 ' 



(5) 



(6) 



7-1 T 

where p is the gas pressure and I is identity matrix. 

The source term on the right hand side of equation (jlj 
is the net gain or loss of energy due to radiative heating 
and cooling and is given (for a static medium) by 



Qr 



ad 



47T 



xT (J* ~ S v ) dv. 



(7) 



This is an integral over frequency v of the difference be- 
tween mean intensity J v and the total source function 
S v , weighted by the total opacity 3 We do not attempt 
to add the corresponding radiation source term to the 
momentum equation. This limits us to applications in 
which radiation pressure is at most a modest fraction of 
gas pressure. An integrator for the coupled MHD and ra- 
diation moment equati ons based on the one-di mensional 
algorithms discussed in ISekora fc Stona (|2010t ) has been 
implemented in Athena and extended to multidimensions 
by JSD12. These more advanced techniques are needed 
to handle the stiff source terms and modified dynamics 
in radiation pressure dominated flows. 

In order to compute the energy source term due to 
radiation, the MHD equations must be supplemented by 
the time-independent equation for RT 

n ■ VI„ = xl ot (Sv -I v ), (8) 
where I v is the specific intensity for an angle defined 
by the unit vector h. In this work, we consider opaci- 
ties due to scattering \ S v an d true absorption x^ bs , with 
xi ot — Xv hs + X S v- It is convenient to define the pho- 
ton destruction probability e„ = x^ bs /xl ot - The source 
function is then given by 

S v = e v B v + (l-e v )J v , (9) 

where B u is thermal source function. The mean intensity 
J„ is the "zeroth" moment, or average, of /„ over solid 
angle 

J v = J Iv(h)dQ. (10) 

When absorption dominates e v — > 1 and S v — > B Vl but 
when scattering dominates e„ — > and S v — > J v - Note 
that this expression assumes that scattering is isotropic. 
Although this is not strictly true for many scattering 
processes (e.g. electron scattering), it will generally be a 
good approximation for problems of interest. 

In addition to J v we will also use H„ and K„, the first 
and second moments, respectively. Their components are 
given by 

Ht = — ( I v {n)pidtt, (11) 



K 1J = ■ 



-in 
1 

4-7T 



Iv(n)pipjd£l, 



(12) 



3 Note that x^ ot has units of [cm 1 ]. Throughout this work 
we will use \ f° r quantities with these dimensions and k = \l P 
for quantities with dimensions of [cm 2 /g], but will refer to these 
interchangeably as opacities. 
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where dfl is the differential of solid angle, and /Ltj = n • 
Xi. These moments are related to the radiation energy 
density E Ta &, radiation flux F ra( j, and radiation pressure 
Prad via the standard definitions 
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/ J v dv, 
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(13) 


F ra d 


= 47T 


/ H u dp, 
'0 

POO 

/ K v du. 

Jo 


(14) 


Prad 


47T 


(15) 



Integration of equation ([8]) over solid angle yields 



- V • F rad = 4tt 



dv. 



(16) 



and provides an alternative (differential) form for the ra- 
diation source term in equation (j4}. The differential form 
tends to perform better in regions where optical depths 
across a gridzone are large, while the integral form is 
preferable in regions of low optical depth. Hence, we will 
use both expressions, as discussed in section |U 

We have not been forced to make distinctions between 
the Eulerian and comoving frame for radiation vari- 
ables as we have dropped all velocity dependent terms 
in equations $7}, J5J), and (ITH1) . We neglect these terms 
because they are negligible for the tests considered in 
this paper. However, we anticipate solving problems 
where the velocity dependent terms may be important 
and can implement terms that are first order in v/c in 
our RT solver, where necessary. For consistency with 
the VET solver (JSD12), we will adopt the mix frame 
approach where /„, its moments, is, and h are Eu- 
lerian frame variables, while opacities and emissivitics 
are defined in the comoving frame. Derivations of the 
mixed f rame equations can be fou n d inlMihalas fc Klein! 
(fl982l,. iMihalas fc Mihalasl (Il98l . lL owrie et all (I1999L 
and lHubenv fc Burrows! (|2007t) . 

Since we neglect the time derivative of I v and terms 
that arc first order in v/c in equation ([5)1. our method 
is only formally reliable in the static diffusion and free 
streaming-limits. Specifically, the timescale for fluid flow 
if ~ L/v across a characteristic length scale L in the 
simulation domain must be longer than the time it takes 
for radiation to diffuse i dif ~ L 2 \ t ot /c or free-stream 
if s ~ L/c across the domain (see e.g. IMihalas fe Mihalasl 
Il984f) ) . This is sufficient for the test problems considered 
here and should be adequate for many of the problems 
of primary interest to us. When necessary, we can retain 
terms first order in v/c in equations (0 and (jHJ and the 
code will be formally accurate in the dynamic diffusion 
limit (if < idif) as well. 

Throughout this work B v is assumed to correspond 
to the Planck function and is a function only of v and 
gas temperature T. We assume an ideal equation of 
state with p = pRT and gas thermal energy density 
E gas — 2?/ (7 — 1)- Here R is the gas constant and 7 
is the adiabatic index. The adiabatic sound speed is 

a = V IV I P- 

The methods for solving the MHD equations without 
the radiation sourc e term are described in detail in pre- 
vious publications dGardiner fc Sto"nel 120081: IStone et"al"1 
120011 iStone fc Gardinerll2'ool " and are unchanged by the 



solution of radiation transfer. The computation of RT is 
described in Section [3] and the interface of the RT solver 
and MHD integrator is described in Section 2J The se- 
quence for a single timestep can summarized as follows: 

1) Using the hydrodynamic variables (typically T and 
p) from the previous timestep as inputs, we compute x£/ 0t > 
e v , and B v , or each frequency in each grid zone. 

2) We solve Equation ((5J using the methods described 
in Section |3l yielding S v and J v everywhere in the do- 
main. 

3) Using S v and J v (or H v ), we compute the radiation 
source term Q la _d and update equation ((4|) as described 
in Section |H 

4) We advance the MHD variables using the standard 
Athena integrators. 

3. SOLUTION OF RADIATION TRANSFER 

An extensive literature on the solution of RT for 
astrophysical problems in multidimensions exists and 
there are nume rous monographs and rey i ew articles on 
the topic (e.g . IMihalas fc Mihalasl 11984 ICastorl 12004 
iCarlssonl 120081 ). With this literature to draw from, we 
have largely adopted a strategy of implementing existing, 
well-developed algorithms. Since there are many differ- 
ent methods with different strengths and weaknesses, the 
major challenge is finding a method which best suits our 
particular needs. Our most salient constraints include: 

1) The method needs to be amenable to domain de- 
composition since this is the primary algorithm for par- 
allelizing the solution of the MHD equations in Athena. 

2) The method must be able to handle the explicit 
dependence of the source term on J„ in equation Q for 
problems in which scattering is important (i.e. we must 
be able solve non-LTE problems). 

3) The method needs to be able to handle (shearing) 
periodic boundary conditions. 

4) The method must be robust and capable of handling 
discontinuities in temperature and density which arise 
when shocks are present in the flow. 

5) Ideally, the method should be efficient enough that 
for simple problems (e.g. LTE with grey or mean opaci- 
ties), neither the memory constraints nor the total com- 
putational time is dominated by the solution of RT. 

With these considerations in mind, we have 
implemented a short -cha racteristics base d solver 
(IMihalas et all 119781: lOlson fc Kunaszl Il987t 
iKunasz fc Auerl Il988f h In this method the specific 
intensity is discretized on a set of rays at each cell center 
in the simulation domain. Equation ([5} is integrated 
along each ray using initial intensities interpolated 
from neighboring grid zones. Since only neighboring 
grid zones are used for this integration, the total 
computational cost (per iteration) scales linearly with 
the number of gridzones in the domain. This is also 
simple to parallelize with domain decomposition as only 
information from cells on the faces of the neighboring 
sub-domains need to be passed. 

Thi s is in contra st to a long characteristics method 
(e.g. lFeautrierlll964l ). which would integrate the RT equa- 
tion along each ray through all gridzones intersected by 
the ray until the edges of the simulation domain are 
reached. Such a method is generally more computation- 
ally expensive since computation of the specific intensity 
in each gridzone typically requires integrating equation 
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iJSJ through ~ iV 1 / 3 gridzones (where N is the total num- 
ber of gridzone in the domain). It is also more cumber- 
some to use with dom ain decomposition (see, however, 
iHeinemann et aT1l2006l ) since it may require the passing 
of larger blocks of data, including information from non- 
neighboring subdomains. 

Although short characteristic methods are computa- 
tionally more expedient, they suffer from greater nu- 
merical diffusion due to the interpolation that is re- 
quired to compute the intensity in neighboring gridzones 
(|Kunasz fc Auerll"l988l ). For problems where a few grid- 
zones (or point sources) dominate the total emissivity, a 
short characteristics solver may require very high angu- 
lar resolution to accurately resolve the radiation field far 
from the dominant source. If the angular resolution is 
too low, anomalous structure (e.g. spokes) in the heat- 
ing and cooling rates will emanate from the dominant 
sources (see e.g. iFinlator et all 12009). (In this case the 
numerical diffusion introduced by interpolation can be 
beneficial.) Instead, the emission from point sources is 
better ha ndled by suitably designed long characteristics 
methods (lAbel fc Wandeltll2002l; iKrumholz et aTll2007bl 
although see also lRiikhorst et al.ll2006l ). For the applica- 
tions of interest in this work (e.g. accretion flows), the 
diffuse radiation field dominates. Moreover, even when 
point sources are present, the diffuse radiation field due 
to scattering or re-emission (e.g. HII regions) cannot 
generally be ignored, and therefore we anticipate such 
problems may be accommodated in the future by a hy- 
brid scheme which uses short characteristics for the dif- 
fuse emission, and long characteristics for bright point 
sources. 

Non-LTE problems are handled via iteration. For 
each time step the formal solution of the whole do- 
main is repeated, updating J v and S v during each 
iteration, until some formal convergence criterion is 
met. As discussed below, we implement an accel- 
erated (or approximate) lambda iteration (hereafter 
A LI) algorithm based on the Gau s s-Seid el method 
of iTruiillo Bueno &; Fabiani Bcndicho (1995) (hereafter 
TF95). The TF95 method is efficient for solving non- 
LTE problems because it significantly increases the con- 
vergence rate without significantly increasing the com- 
putational cost (or memory footprint) per iteration. 

Iteration is also used in LTE problems to handle 
boundary conditions at the interface of subdomains and 
for physical periodic boundary conditions at domain 
edges. On each iteration the incoming intensity from the 
neighboring subdomain is fixed from the previous itera- 
tion (or timestep for the first iteration). The outgoing 
intensity, which corresponds to the incoming intensity 
in the neighboring subdomain, is then updated and the 
formal solution is iterated to convergence. For LTE prob- 
lems, this is not the most efficient me thod for handling 
the s ubdomain boundaries (see e.g. IHeinemann et all 
2006). For the moment, we are primarily interested in 
non-LTE problems where iteration is required regardless. 
We generally find fairly rapid convergence (requiring only 
a few iterations) for most of our LTE test problems when 
iterations is used, so this is not a significant limitation. 

In many respects our s hort-characterist ic s RT solver 
is similar to tho se of Ivan Noort et al.1 <|2002f ) and 
lHavek et al.l ((2010) in that both implement ALI to han- 
dle deviations from LTE and both utilize domain de- 



composition for parallelization. lHavek et al.l (|2010T) used 
their code to solve the RT equation including scatter- 
ing, in MHD simulations of stellar atmospheres on three- 
dimensional Cartesian grids. Hence, the effectiveness 
of several key aspects of our module have already been 
demonstrated in a sophisticated MHD code and applied 
to realistic astrophysical applications. 

3.1. Frequency Discretization 

The scheme we have implemented allows for the com- 
putation of frequency dependent, grey, or monochro- 
matic RT. Radiation variables (moments and specific in- 
tensities) and radiative properties of the fluid such as 
the opacities, thermal source function, and photon de- 
struction parameter are tabulated on a grid of ri{ > 1 
discrete frequencies or frequency groups. For flexibility, 
the functional form of opacities and emissivities can be 
specified via user-defined functions. In general, the com- 
putational cost and memory footprint of problems scale 
linearly with Uf. 

These frequency bins can simply be discrete frequen- 
cies when RT is used to generate diagnostic outputs such 
as images an d spectra. Group mean opacities and emis- 
sivities (e.g. iMihalas fc Mihalasl Il98l ISkartlienl [2000) 
and corresponding quadrature weights must be specified 
when the RT solver is used to compute the radiation 
source terms or VET. In the simplest case, rif — 1 and an 
appropriate frequency integrated mean opacity is speci- 
fied. 

Unless otherwise noted, we will drop subscripts de- 
noting the frequency dependence of radiation variables 
and only describe the monochromatic problem hereafter. 
For the problems under present consideration, there is 
no explicit coupling of the specific intensity at different 
frequencies so the frequency dependent calculation is a 
trivial generalization of the monochromatic problem. 

3.2. Angular Discretization 

We discretize the specific intensity on both angular 
and spatial grids. For one-dimensional problems, the dis- 
cretization is chosen so that polar angles correspond to 
the abscissas for Gaussian quadrature. In multidimen- 
sions, discretization of the angles proceeds according to 
the al gorithm described in Appendix B of iBruls et al.1 
(1999), which is b ased on the prin ciples of type A quadra- 
ture described in ICarlsor] ([1963). 

This method attempts to distribute the rays as evenly 
as possible over the unit sphere, subject to the constraint 
that each octant of the unit sphere is discretized identi- 
cally. Hence the angle discretization is invariant for 90° 
rotations about the coordinate axes. This is desirable be- 
cause Athena is designed to be a general purpose code, 
and there is often no preferred direction with which to 
align the angular grid (as in some atmosphere calcula- 
tions). Without this constraint, the result would gener- 
ally depend on the orientation of boundary and initial 
conditions relative to the coordinate axes. 

The user specifies the number of polar angles n^, and 
the code generates an array of n a rays covering the unit 
sphere. In one dimension, this corresponds to n a = 
rays because of axisymmetry. For multidimensional do- 
mains n a = n^(n^ + 2) rays. However, in two-dimension 
only half of these are unique due to the implied invari- 
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Fig. 1. — Schematic of the RT solution for an individual grid- 
zone whose cell center correspond to vertex E in a two-dimensional 
radiation grid. In this case S°, and x° are known at E, and 1° 
is to be computed. 1° is computed along a each ray using equa- 
tion H20II . Since quantities I~ , S^, and x do not correspond to 
vertices of the grid, they must be computed via interpolation from 
neighboring grid zones. The red ray intersects rows between ver- 
tices A and B (upwind) or H and I (downwind). Hence the values 
of S, X an d I at these vertices are used to interpolate S 1 , x , 
and I~ . Vertices C and G are also used with quadratic interpola- 
tion. A similar case holds for the blue ray, but with interpolation 
performed on columns A-D-G and C-F-I. The open and closed cir- 
cles denote vertices which are used for the interpolation of I~ . The 
dashed curve is an extension of the blue ray which intersects A-B-C 
row. 



ance of physical quantities in the third dimension and 
n a = n. M (n M + 2)/2. 

Setting — 2 in a one-dimensional calculation is 
analogous to invoking the two-stream approximation, in 
which the radiation field of each hemisphere is approxi- 
mated by transfer along a single ray. This assumption is 
commonly used to derive analytic solutions, and allows 
the ratio of H/J to vary but keeps the ratio of K/ J fixed 
at 1/3, consistent with the Eddington approximation. In 
two (three) dimensional calculations, choosing = 2 
approximates each quadrant (octant) with a single ray 
and also yields = 1/3(5™ J. The algorithm is wel l- 
defined and unique only for n M < 12 ()Bruls et alj fl999). 
corresponding to n a = 168 (84 in two dimensions). This 
should not be prohibitive for the problems of interest. 

For each ray fik, we compute a vector of direction 
cosines (nok, /J-ik, M2fe) with fj,^ = hk ■ Xi and quadrature 
weights Wk- Then equations f(TTJ|) - (jT2 |) become 

Tla-l 

j=Y, w ^ ( l7 ) 

k=0 
n a -l 

Hi = ^2 Wkh^ik (18) 

k=0 
n a -l 

Kjj = \] WklkfijkfJ'jk, (19) 

k=0 

where I k = I(fik)- 

3.3. Implementation of the Short- Characteristics 
Algorithm 



The short characteristic method (iMihalas et al.l 119781 : 
lOlson fc Kunaszl 119871; iKunasz fe Auerl 119881 ) has been 
discussed previously by several authors. The basic com- 
putation step for a single gridzone in both LTE and 
non-LTE problems is illustrated in Figure Q] for the two- 
dimensional case. Fluid radiative properties and radia- 
tion variables (e.g. x tot , B, e, Ik, J, S) are defined on a 
radiation grid. The vertices of this grid correspond to the 
cell centers of the MHD domain so that fluid radiative 
properties are computed directly from the cell centered 
MHD fluid variables. Generalizing to three dimensional 
domains is straight-forward. 

At each vertex, the specific intensity at x° is com- 
puted along each ray hk from x^ to x£ . For second-order 
interpolation the intensity is given by 

I°k=Ik ei ~ ATn + *kSk + **<S° + *^, + , (20) 
where , ^F", and Vf^ denote interpolation coefficients 
which depend on the opacities xZ", x° , and xt through 
the optical depth intervals At^ and At^ . 

The form of the interpolation coefficients fyZ, ^>1, and 
Vl/^ depends on the interpolation method used. The stan- 
dard expressions for secon d order interpola t ion ar e listed 
in equations (7a)-(9c) of IKunasz fc Auerl ([1981 . One 
drawback of these expression is that they are subject to 
overshoot where gradients in S^ and Xk are st ee P- For- 
tunately, these cases can b e handled wi th B ezier-type in- 
terpol ation as described in I Auerl (|2003l ) and lHavek et al.l 
(2010). With Bezier-type interpolation schemes, one 
can utilize a control point S% = S° — 0.5AtZ (dS/dr)° 
to determine if overshoots are present in the standard 
second-order expressions. If S^ < mhx(Su , S°) or S% > 
max(SZ , S°) overshoots are present and alternative ex- 
pressions are utilized. Suit able choices f o llow f rom set- 
ting S c k = SZ or SI = S°. lHavek et all (|20ll provide 
the corresponding expressions for vE^ and in their Ap- 
pendix A. Similar methods are also used to compute in - 
tervals At^ (see e.g. equation A. 3 of lHavek et alJ feOlO). 

For one-dimensional problems, and x^" correspond 
to neighboring grid vertices Xi-\ and Xi+\. Hence, 
IZ = 7i_i,fe, which was just computed in the neighboring 
zone while S^, and Xk can be computed directly from 
hydrodynamics variables at sc,*±i. In multidimensional 
problems, x^ and no longer correspond to vertices of 
the radiation grid and variables IZ , Sj^ , and Xk mus t 
be interpolated. We implement and test both first-order 
(linear) and mon otonic second order (q uadratic) inter- 
polation schemes (|Auer fe Paletoulll994D . Both methods 
prevent overshoots and enforce positivity of the inter- 
polants. The choice is particularly relevant for Ik, as 
second-order methods generally produce much less diffu- 
sion of the radiation beam. A drawback of second order 
interpolation is that it places additional constraints on 
the order in which one sweeps through gridzones and the 
stencil used for the evaluation of Ik ■ 

Consider the two rays depicted in Figure [TJ We com- 
pute interpolants S^ and xt using known quantities at 
vertices of the radiation grid. If row A-B-C or column 
A-D-G correspond to ghost (boundary) zones, IZ can be 
computed from the (prescribed) boundary intensities. If 
they are not ghost zones, interpolation can only be per- 
formed on zones in which Ik has already been computed. 
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Fig. 2. — Progression of the sweep through a two-dimensional 
grid (domain or subdomain) when linear interpolation is used. The 
forward sweep (red curve) first progresses parallel to y, computing 
RT only for upward pointing rays (ft ■ y > 0). For each row (fixed 
yA, one first sweeps parallel to x, computing RT along rays with 
ft • x > until reaching the grid edge rrjv, then reverses direction 
and computes along rays with ft ■ x < until reaching the grid edge 
xq . This continues until one reaches gridzone (xq, y^ ) . The back- 
ward sweep (blue curves) inverts the forward sweep, computing RT 
only for downward pointing rays (ft • y < 0). In the Gauss-Seidel 
method, updated values of Si j are incorporated into the back- 
ward sweep, beginning with Sm,n- The three-dimensional case is 
a straightforward generalization. 



If we first sweep along rows of fixed yj (as in Figure [2]), 
Ik has only been computed at vertices A, B, C, and D. 
This means that Ik is known for all vertices used in the 
linear interpolation of Ij~ as well as for quadratic interpo- 
lation (and any higher order interpolation) of rays which 
intersect row A-B-C. 

We refer to rays that intersect the column A-D-G, such 
as the blue one in Figure [IJ as "shallow" rays. Shallow 
rays are a potential problem for quadratic (and higher 
order) interpolation, since Ik at G has not been com- 
puted. When quadratic or higher order interpolation is 
desired, such rays can be handled in a number of ways. 
One possibility is to switch the order of the sweep for 
shallow rays so that it first proceeds in the y direction 
along columns of fixed Xi. In this case Ik for shallow 
rays will be known at vertices A, B, D, and G. The main 
drawback (discussed further in Section EOl below) is that 
one is unable to implement a Gauss-Seidel iteration for 
non-LTE problems. 

One can also construct alternatives by extending the 
stencil beyond vertices A-I. For example, one can extend 
shallow rays until they intersect row A-B-C as shown by 
the dashed curve in Figured] A drawback of this solution 
is that it requires modest additional effort for computing 
\&7T, although this can be alleviated by computing only 
on the first iteration and reus ing it for subsequent iter- 
ations (e.g. lHavek et al.ll2010t ). Alternatively quadratic 
interpolation could be preformed using A, D, and the 



vertex directly below A ()Kunasz fc Auerlfl988l) . 

These two solutions share common drawbacks. For 
parallelization with domain decomposition, only one 
ghost zone is needed per grid zone on a subdomain face, 
when only vertices A-I are used. Extension of shallow 
rays beyond this stencil requires the passing of additional 
data and associated bookkeeping. More philosophically, 
we feel it is desirable to treat all rays as consistently 
as possible. In either of these schemes, RT along some 
rays will be computed using only neighboring grid zones, 
while other rays will not. Our preference is to treat all 
rays on the same footing. 

For this reason, we have decided to switch the order 
of the sweep for shallow rays. Athena is implemented 
so that each sub-grid of the domains has regular spacing 
and therefore gridzones with fixed aspect ratio. This 
means that the distinction between rays that are shallow 
and those that are not is equivalent for each grid zone. 
However, our definition of a shallow ray depends upon 
the direction of the sweep. The blue ray in Figure [T] is 
shallow because we first traverse the grid along rows of 
fixed yj, only moving to yj+i when intensity has been 
computed for all gridzones in the row yj, as depicted in 
Figure [H 

If we reverse the sweep so that we first traverse columns 
of fixed Xi, the blue ray will no longer be shallow, as the 
intensity at G will be computed before it is needed for 
the computation of the intensity at E. In this case the 
red ray is now a shallow ray as the intensity at C will not 
have been computed before it is needed to compute the 
intensity at E. Hence, by varying the sweep direction, 
we can handle all rays and accommodate a quadratic 
interpolation scheme which computes all intensities in a 
gridzone (xi,yj) only using intensities from neighboring 
gridzones (x i ±i,y J ±i). 

3.4. Iterative Methods for Non-LTE Problems 

We now describe how we handle non-LTE problems it- 
eratively. Following common convention we denote the 
angle averaged formal solution of the RT equation (here- 
after, simply the formal solution) in operator notation 
as 

J = AS. (21) 

Here, A is a linear operator representing the (discretized) 
formal solution, and J and S are vectors spanning each 
gridzone in the simulation domain. Using equation ([9]) 
to eliminate J, one obtains an equation for S in terms of 
B 

S = (1 - e)A[5] + eB. (22) 

Since A is a linear operator we can solve for S 

S=[l-(l~e)A]- 1 [eB}. (23) 

If one can invert A a formal solution of the non-LTE 
problem follows from solving (|23j) and obtaining J from 
(|2I[) . However, for three-dimensional problems A is a 
very large matrix and not sparsely populated when sys- 
tems are far from LTE so its direct inversion is imprac- 
tical. Therefore, equation ([22]) is usually solved via iter- 
ation. 

A simple iterative scheme for solving equation (|22|) be- 
gins with an initial guess for the source function S , 
which is then used to compute an improved estimate 
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S n+1 = (1 - e)A[S n ] + eB. However, this method (often 
referred to as Lambda Iteration) has very poor conver- 
gence properti es. For practical pro blems, ALI methods 
Cannon! 1 19731) a re com monly used. iRvbicki fe Hummer! 



19911 ). iHubenvl ([2003h and TF95 provide useful reviews 



of ALI methods and we refer the reader to these works 
for a more in-depth discussion. Here we just summarize 
the basic concepts involved. 

In ALI methods one solves equation (|23|) directly, but 
using an approximate form A* which is easier to invert 
then the full A operator. Since only the approximate 
A* is used, iteration is still necessary. Numerous choices 
for A* have been proposed, but it has been argued that 
simply taking the diagonal elements of the full A ma - 
trix represents a n ear-op timal choice (|01son et al.lll986l ). 
lOlson fe Kunaszl JI987) provide expressions for diago- 
nal elements of A when short characteristics are used. 
In each grid zone the change in the source functions 
AS, = S? +1 - S" can be written as 



AS, = 



(1-e.QJf+e^-fff 
1 - (1 - e i )A ii 



(24) 



where the subscript i enumerates all gridzones in the do- 
main. 

As TF95 discuss, when J™ is exclusively used in equa- 
tion (|2"4"|) . the ALI scheme is equivalent to the Jacobi iter- 
ative method for solving linear systems. TF95 show that 
one can construct a Gauss-Seidel algorithm by incorpo- 
rating the new values of in equation (|2"4l as these 
become available. Here i' < i refers to gridzones where 
J has already been updated. The complexity of devis- 
ing a Gauss-Seidel algorithm for RT comes from the fact 
that the computation of specific intensity I^k for some 
subset of the rays hk need to be computed using S™ r<i 
rather than (i.e. old rather than new values of the 
source function). Therefore the contribution from these 
particular rays to * must be corrected as the updated 

values S!pw become available. 

TF95 give a detailed description of how to implement 
such an algorithm on a one-dimensional domain. The al- 
gorithm requires storing a modest amount of data in each 
gridzone, but very little additional computation. The 
convergence rate is improved by a factor of two, so prob- 
lems requiring several iterations gain nearly a factor of 
two decrease in computational effort for only a minor 
increase in code complexity. 

When linear interpolation is used, the generaliza- 
tion of their one-dimensional method to two and three- 
dimensional domains is straight-forward. The two- 
dimensional sweep proceeds as depicted in Figure [5J The 
vertices in the radiation grid correspond to cell centers 
(xi,yj). The sweep generally proceeds with i as the 
more rapidly varying index. Consider a domain with 
N x = N y = N for simplicity. In each gridzone (xi , j/j ) , 
we first compute the intensity Ik for all upward directed 
rays (hk ■ y > 0) in the forward sweep and then for all 
downward directed rays (hk-y < 0) on the reverse sweep. 

On the reverse sweep, the upper right gridzone 
(x./v, j/ at) is the first in which the computation of all new 
intensities is completed. At this point J^ + ^ is com- 
pletely specified and we compute S^^. From here on, 
all subsequent RT computations use S^^ rather than 



These must also 
N N and weights 



&n n- However, this alone is not sufficient to make 
it a Gauss-Seidel scheme, because the contributions to 
Jff—i n> ^nn—v an d JfftSi ff—i lrom upward directed 
rays on the forward sweep used S 1 ^ N . 

be updated using ASn,n — S^~ N — S' 
which were saved on the forward sweep. We also up- 
date the outgoing intensities I^ +l (since they were also 
computed using S 7 ^ N ) as they correspond to the incom- 
ing intensities in neighboring gridzones. Since the corre- 
sponding weights have already been computed as part of 
the forward sweep, the additional computational cost is 
very modest. 

Following the discussion in Section 13.31 we note that 
feasibility of performing a Gauss-Seidel iteration with 
quadratic interpolation is dependent on the way shal- 
low rays are handled. Reorienting the sweep for shal- 
low rays so that j is more rapidly varying index, but 
keeping i as the rapidly varying index for remaining rays 
does not allow for an efficient Gauss-Seidel scheme be- 
cause some of the necessary J n+1 (and therefore S n+1 ) 
are not available when the backward sweep begins. 4 In 
the light of this issue, we have implemented Gauss-Seidel 
routines only with linear interpolation. For problems 
where quadratic interpolation is preferable, we default 
to the Jacobi method (i.e. standard ALI). 

We continue the iteration until some convergence cri- 
terion is met. Consistent with previous work, we stop it- 
erating when the maximum relative change in the source 
function over the whole domain is less than some pre- 
scribed threshold S c 



AS, 
S, 



(25) 



For LTE problems that use iteration to handle boundary 
conditions, S does not change from one iteration to the 
next and we replace S% with Ji in equation (|25[) . 

The choice of 6 C is clearly an important input to the 
method, but there is no firmly established criterion and 
the optimal choice depends on a number of considera- 
tions that may be problem dependent. Since the compu- 
tational cost of the method generally scales linearly with 
the number of iterations performed and a lower threshold 
leads to more iterations, there is a tradeoff between ac- 
curacy and computational expediency. With the excep- 
tion of the uniform temperature non-LTE atmosphere, 
the tests presented in section [5] were performed using 
5 C = 10~ 5 . Increasing S c to 10~ 4 had a negligible impact 
on the linear wave tests. 

Our expectations based on the tests we have performed 
so far are that for the problems of primary interest to us 
(e.g. shearing box simulations of accretion disks) S c ~ 
10~ 5 — 10~ 3 will be suffic ient, consistent w ith studies 
using similar methods (e.g. lHavek et al.ll2010t ). However, 
we emphasize that the appropriate choice will be problem 
dependent and must be assessed on a case-by-case basis. 
We view the choice of 5 C in roughly the same terms as 
we view the choice of grid resolution. One can adopt a 
threshold based on previous results and experience, but 
ultimately one needs to compute the problem using a 

4 For two-dimensional domains one can devise an efficient Gauss- 
Seidel algorithm that sweeps diagonally through the grid, but this 
implementation does not generalize to three dimensions. 



range of S c and choose a sufficiently small value such 
that the results are insensitive to the choice. 

3.5. Boundary Conditions and Parallelization 

Boundary conditions and domain decomposition in 
Athena are both implemented for MHD via the use of 
ghost zones, and we implement RT boundary conditions 
in an analogous way. The solver computes RT in grid- 
zones on a boundary (domain or subdomain) in the same 
way as an interior gridzone, but using the intensities and 
source functions from the ghost zones to compute the 
relevant integration weights and interpolants. The in- 
tensities and source functions in the ghost zones are de- 
termined according to prescribed boundary conditions. 

In general, boundary conditions for the MHD integra- 
tor will not translate directly to boundary conditions for 
the RT solver. Different problems with the same MHD 
boundary conditions may require different boundaries for 
the radiation field. Hence separate boundary conditions 
must be prescribed when using the RT solver. For the 
code test problems presented in section [5j we have im- 
plemented two types of boundary conditions specifying 
either fixed incident intensity or periodic intensities on 
the boundaries. Other boundary conditions can be spec- 
ified via user defined functions. 

Athena runs on parallel machines using domain decom- 
position implemented through MPI calls. The MHD inte- 
grator passes all conserved variables and passive scalars 
from faces of neighboring subdomains to ghost zones. 
The MHD integrator requires either four or five ghost 
zones for each gridzone on the subdomain face. The 
RT solver operates analogously, passing intensities and 
source functions, but only requires one ghost zone for 
each gridzone on a subdomain face. 

The main differences between the RT solver and the 
MHD integrator are the frequency and quantity of data 
that must be passed. For each frequency bin in every 
ghost zone I must be passed for all n a rays with quadratic 
interpolation or, alternatively, n a /2 incoming rays with 
linear interpolation. For non-LTE problems S and H 
must also be passed. Hence for quadratic interpolation, 
the code passes a total of n/(n a + 1 + ndi m ) floating point 
variables per face gridzone per iteration, where ridim is 
the number of dimensions in the domain. In contrast, 
the MHD integrator typically passes ~ 50 floating point 
variables per face gridzone per timestep. For problems 
where n a and n t are small and few iterations are required 
(e.g. an LTE grey problem), the volume of RT data is 
therefore comparable to and may even be less than the 
amount of data passed by the MHD integrator. 

We note that the use of iteration to handle subdo- 
main boundary conditions may lead to some dependence 
on the number of subdomains that are used. We have 
considered the sensitivity of our results to this issue by 
performing most of the tests described in section [5] both 
with and without domain decomposition. In practice, 
the converged mean intensities do not differ (relative 
to the non decomposed domain) by more than ~ S C J. 
The sensitivity is highest for problems where the optical 
depth across an individual subdomain is of order unity or 
smaller, Problems with optically thick subdomains gen- 
erally lead to smaller discrepancies. Since we already 
choose our convergence criterion to be at a level that 
minimizes the impact on our results, this sensitivity to 



the domain decomposition should not lead to significant 
errors. 

4. INTERFACE OF THE RADIATIVE TRANSFER SOLVER 
TO THE MHD INTEGRATOR 

There are two regimes in which the effect of radia- 
tion on the MHD is important. The first is when the 
radiation field is a significant contribution to both the 
energy and momentum fluxes in the flow. In this regime, 
the radiation source terms in the MHD equations can be 
very stiff, and the equations contain wave modes which 
propagate at close to the speed of light. Both of these 
properties require significant modification to the under- 
lying MHD integrators in order to enable accurate and 
stable integration. In JSD12 we describe a method for 
this regime based on an extension of the modified Go- 
dunov method of SS10 to multidimensions, with a VET 
(defined as f = Prad/^rad) computed from a formal so- 
lution of the RT equation using the module described in 
this paper. At each time step, the RT solver computes 
the radiation field as described in Section |31 evaluating K 
and J via equations (fl7|) and ([T9l) . We the compute the 
VET using f = K/J as described in section 3.4 of JSD12. 

The second regime is when the radiation pressure can 
be ignored, and the effect of radiation is only through the 
heating and cooling source terms in the energy equations. 
In principle, the modified Godunov method adopted in 
JSDI2 would be an attractive approach for handling the 
stiff energy source term that can arise in this regime as 
well. However, the modified Godunov method requires 
that one compute the gradient of radiation source terms 
on the plane of primitive variables. This in turn requires 
analytic expression for the radiation sources in terms of 
the fluid variables. Hence, it is generally not a viable 
method for problems where the radiation properties are 
complicated functions of frequency and fluid variables, 
as may be the case with bound-free and bound-bound 
atomic opacities or Compton scattering. 

These limitations motivate us to implement an alter- 
native method to directly compute the radiation source 
term in the fluid energy equation (j4} and couple it to 
the standard MHD integrators. When operating in this 
mode, we perform the formal solution at the beginning 
of each timestep. We first compute fluid radiation prop- 
erties in each gridzone Xj of the domain. This includes 
the variables xT*' ^ an ^ e ^ which are computed via 
user-defined functions of the conserved MHD variables 
and passive scalars from the previous timesteps. We use 
these, along with Ji from the previous time-step, to ini- 
tialize Si. Once the formal solution is completed, we 
account for the source function on the right hand side of 
equation (0| via an operator split update of E. We first 
compute the radiative source function in each zone and 
then update the total energy 

AEi = 6t(Q I3d )i. (26) 

The standard MHD integration algorithm then proceeds 
using this "new" value for E{. 

We compute Q r ad in one of two ways, depending on the 
characteristic optical depth. We either use the integral 
form 

Q'f = *mc? t (Ji ~ St) = ^ X f a (Ji - B t ), (27) 
or the differential form 

Qf f = -4ttV • H i5 (28) 
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Previous work (|Bruls et al.l ll999. and references therein) 
has demonstrated that the integral form is inaccurate 
when the optical depth per gridzone is large. In this case 
Ji — Bi <C Bi while xf° s is large so round-off errors can 
be greatly amplified. The int egral form, howev er, is more 
accurate when xf bs Axi < 1 (jBruls et al.lll999l ). 

Therefore, we have designed our RT solver to com- 
pute either form of Q ra d, depending on the regime of the 
computation. In most applications of interest, there is a 
transition from optically thick to optically thin regions, 
so we must specify a criterion for switching between the 
differential and integral forms in the same domain. For 
the test problems considered here, we find a simple switch 



(Qrad)z 



Qtat if x tot ^ < X 

Qf lf otherwise 



to be sufficient. This has the advantage that it is a purely 
local criterion. Using a method which more smoothly in- 
terpo lates between the two regimes (see e.g. lHavek et al.l 
|2010| ) did not improve performance in a measurable way, 
but may be preferable for more sophisticated applica- 
tions. 

Due to the explicit update, we must take care in choos- 
ing a time step. In the absence of radiation the MHD 
integrator chooses a time step Stc based on the CFL con- 
straint. In principle, this time step can be much larger 
than the radiative cooling time, which could lead to obvi- 
ous errors, such as the energy density becoming negative. 
As we elaborate upon in section lS^l one can derive a gen- 
eralized CFL condition for a radiating fluid based on the 
need to resolve the damping time for a non-equilibrium 
radiation diffusion mode. This time scale <5i r d is gener- 
ally most restrictive when the optical depth per gridzone 
X tot Ax ~ 1, in which case 



St 



rd 



E gas a 
~? dt<i 

Aad C 



(29) 



assuming the adiabatic sound speed a (rather than the 
Alfven speed) sets the CFL condition. This generalized 
CFL constraint can be quite stringent, requiring short 
time steps and increasing computational costs if either 
a <C c or i? ra d > E gas . Hence, many problems will re- 
quire the use of the VET method described in JSD12, 
which uses timesteps determined by the standard (non- 
radiative) CFL constraint. In practice, we are almost 
always limited to problems with E giLS < E rat \, so we do 
not attempt to include the radiation momentum source 
term in equation ^ as it is generally small for problems 
that are computationally feasible with operator splitting. 

The algorithm described above will, in general, only 
be first order convergent. Note that we could construct 
a second-o rder scheme when usin g Athena's VL+CT 
integrator ([Stone fc Gardiner] I2009D . by performing the 
operator split update before the corrector step in the 
predictor-corrector scheme. However, some of the ad- 
vantages of the second order convergence will be lost due 
to the increased diffus ivity of the VL+CT re lative to 
the CTU+CT scheme (IGardiner fc Stondl200ll . Hence, 
we have not yet pursued the possibility although it may 
prove to be a useful avenue for future work. 

5. TESTS 

Our test problems fall into roughly two categories: 
stand-alone tests of the RT solver on fixed domains and 
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Fig. 3. — Convergence of the ALI methods on a highly non- 
LTE uniform atmosphere with e = 10 -6 . The curves show the 
maximum relative difference between the numer ical ly computed S 
and the analytic solutions implied by equation (1306 for the Jacobi 
(solid) and Gauss-Seidel (dotted) methods. We compute the nu- 
merical solutions using a one-dimensional domain with 9 gridzones 
per decade in optical depth. 



tests of the coupled MHD integrator and RT solver in 
fully time dependent calculations. The former are par- 
ticularly useful for evaluating the RT solvers performance 
on multidimensional and non-LTE problems. For the lat- 
ter, we focus primarily on simpler LTE problems, so we 
can compare the simulations result directly to precise 
analytic or semi-analytic solutions. 

Further tests of the RT solver as part of the VET 
method are presented in JSD12. 

5.1. Uniform Temperature Non-LTE Atmosphere 

We begin by solving the monochromatic RT problem in 
a uniform temperature, one-dimensional scattering dom- 
inated atmosphere. This test is particularly useful for 
evaluating the RT solver's performance on a non-LTE 
systems and evaluating the convergence properties of Ja- 
cobi and Gauss-Seidel iterative schemes. We adopt the 
two-stream approximation for the RT solver so we can 
compare directly with analytic solutions based on the 
Eddington approximation. Since we assume a uniform 
opacity k and temperature T, the analytic solution is 
only a function of optical depth dr = xdz, the thermal 
source function B and photon destruction probability e. 
With these assumptions the mean intensity J is given by 



/3€T 



J = B 1 



1 + v 7 ^ 



(30) 



We assume that x oc p and p increases exponentially (but 
keep e constant) with distance from the upper boundary, 
which has no incoming intensity. This provides an expo- 
nential variation in r which is well-suited for resolving 
the transition from LTE to non-LTE within the atmo- 
sphere. 

Figure shows the convergence of the true error of the 
numerically derived solutions. This is evaluated as the 
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Fig. 4. — Comparison of numerical (solid) and analytic (dashed) 
solutions of the source function monochromatic, uniform non-LTE 
atmospheres as a function of optical depth. Each set of curves cor- 
responds to a different photon destruction parameter, running from 
e = 10 — 2 to 10 — 10 from top to bottom. We compute the numerical 
solutions simulations using cubic three-dimensional domains with 
64 3 gridzones, distributed over 64 MPI subdomains. The optical 
depth variation is aligned with the z-axis of the simulation domain 
and the solution is uniform in the horizontal directions. 



maximum relative difference \AS\/S, with AS the differ- 
ence of the numerically derived S from the analytic so- 
lution. We first consider a one-dimensional domain with 
e = 10 -6 , as this gives a highly non-LTE atmosphere 
and facilitates direct comparison with Figure 3 in TF95. 
We initialize the radiation field to be in LTE everywhere 
(J = 13). We consider two different iterative schemes: 
Jacobi and Gauss-Seidel. As expected, the convergence 
rate of the Gauss-Seidel methods is nearly a factor of 
two better than Jacobi. We assume nine gridzones per 
decade in t to match TF95 and our convergence rates 
agree reasonably well with those shown their Figure 3. 

We have also implemented the successive over- 
relaxation (SOR) method of TF95, and find rapid con- 
vergence, consistent with that shown in Figure 3 of 
TF95. We have tested SOR on both one-dimensional and 
two-dimensional domains and find that it is an effective 
method as long as all boundary intensities are fixed dur- 
ing iteration. However, if the intensities on one of the 
boundaries vary from one iteration step to the next, the 
method is generally not stable. For example, instabil- 
ity occurred when we used periodic boundary conditions 
or when we employed subdomain decomposition. Since 
most of our primary science goals involve problems that 
require the use of periodic boundary conditions or do- 
main decomposition, we do not consider SOR a generally 
viable method for our work. Nevertheless, it may be an 
effective method for a modest sized problem that can be 
run serially with fixed boundary intensities. 

We next consider the same test problem, but use a 
cubic three-dimensional domain with = 2. We align 
the variation of density with the z axis of the domain 
and use periodic boundaries in the horizontal directions. 
Figure HI shows a comparison of the numerical and ana- 
lytical solutions for various choices of e. The agreement 
between the numeric and analytic solutions is quite good 
overall, but tends to be poorest at low optical depths. 
For fixed resolution, the discrepancies with the analytic 
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Fig. 5. — Comparison of Athena (diamonds) and Feautrier 
(solid) emission spectra from the upper boundary of a one di- 
mensional atmosphere. Both calculations assume isotropic electron 
scattering and free-free (Bremsstrahlung) absorption and emission 
for a completely ionized H plasma. The intensities are computed 
using the same angular grid corresponding to abscissas of a 16 
point Gauss-Legendre quadrature of the interval (1,1). The plot- 
ted intensities (from top to bottom) correspond to cosi =0.10, 
0.28, 0.46, 0.62, 0.76, and 0.99. The atmospheres have constant 
temperature (10 6 K) and density which varies exponentially with 
distance, rising from 10 — 6 g/cm 2 at the upper (surface) boundary 
to 10 — 4 g/cm 2 and lower boundary. For comparison, we plot the 
corresponding blackbody at 10 6 K as a dashed curve. 



solution tends to increase as e decreases and the domain 
deviates more strongly from LTE. The accuracy of the 
numerical solution improves with increasing resolution, 
but the number of iterations needed for convergence in- 
creases roughly linearly with resolution. The number 
of iterations required for convergence also increases as e 
decreases. Hence, greater deviations from LTE require 
a greater number of iterations for convergence, as one 
would expect. 

Although some RT problems do require explicit fre- 
quency coupling (e.g. Compton scattering, partial re- 
distribution), many problems can be treated in the ap- 
proximation that frequencies are not explicitly coupled. 
Multifrequency problems are then just a series of sin- 
gle frequency calculations and, hence, a straightforward 
generalization of the monochromatic problem. Figure [5] 
shows the intensity spectrum from a multifrequency cal- 
culation done with Athena for a uniform temperature at- 
mosphere. We again assume p varies exponentially with 
distance, rising from 10 -6 g/cm 2 at the upper boundary 
to 10~ 4 g/cm 2 and lower boundary. The results plotted 
here are for N x — 256. 

Incoming intensity at the upper boundary is assumed 
to be zero and I v = B v at the lower boundary. We in- 
clude isotropic electron scattering opacity and free-free 
(Bremsstrahlung) emission and absorption. The electron 
scattering is modeled as isotropic and the cross-section is 
the Thomson cross-section. For simplicity free-free pro- 
cesses are computed assuming a Gaunt factor of unity. 
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Fig. 6. — Mean intensity J due to incident beams propagating through a rarefied, two-dimensional domain. The source function and 
opacity are zero everywhere and the horizontal boundaries are periodic. The boundary condition at the top and bottom of the domain arc 
zero incident intensity except for two gridzones at the base of the domain. In each of these gridzones / is non-zero for a single ray. The 
two rays make the same angle with the y-axis, but are oppositely directed in x, with h ■ y = 0.174 and h ■ x = ±0.628. The left and right 
panels show results from computations with linear interpolation and monotonic quadratic interpolation, respectively. 



The plasma is assumed to be completely ionized Hydro- 
gen. Hence, e„, B V1 and Xv are au functions of frequency. 
However, for an individual frequency the calculations are 
very similar to those described above. The only differ- 
ence is that e v is now a function of depth as well, due 
to the different dependence of scattering and absorption 
opacity on p. 

Figure [5] also shows the results of a Feautrier calcu- 
lation for the same atmosphere using the same angular 
grid. The two calculations generally agree quite well, 
although there is a tendency for the Athena solver to 
give slightly higher intensities for frequencies where the 
spectrum deviates from blackbody. The discrepancy be- 
tween the results is a function of spatial resolution with 
agreement between the two codes improves as the N x is 
increased in the Athena calculation. Calculations on two 
and three dimensional domains (but with density varying 
only in one dimension) yield similar results. 

5.2. Beam Tests in Two Dimensions 

We now consider the propagation of crossing beams of 
radiation, incident on the boundary of a rarefied (B = 0, 
X = 0), periodic domain. This test is particularly useful 
for evaluating the amount of diffusion associated with the 
interpolation schemes for the specific intensity. It is also 
useful testing the performance of periodic and subdomain 
boundary conditions. 

The results for a two-dimensional domain with peri- 
odic boundary condition in the horizontal direction are 
shown in Figure The figure compares a computation 
with linear interpolation to one with quadratic mono- 
tonic interpolation. Our implementation of these meth- 
ods is described in Section 13.31 It is clear from Figure 
[6] that linear interpolation leads to substantially greater 
diffusion of the radiation beam. 

Depending upon the application, the additional dif- 
fusion in the linear interpolation scheme can be cither 
advantageous or problematic. On one hand, a less diffu- 
sive scheme allows one to model important effects, such 
as shadowing by optically thick material, with greater 
fidelity. Indeed, the ability to more accurately capture 



such effects is an important motivation for using RT in- 
stead of more ad hoc closure prescriptions, such as FLD. 

However, computational expedience limits the angular 
resolution we can achieve. When only a modest number 
of rays are used with a less diffusive scheme, fan-shaped 
"spokes" can appear in the mean intensities and Edding- 
ton factors, if the emission in a small number of grid 
zones significantly exceed that of surrounding zones. In- 
deed, our short characteristics based method is not well 
suited to problems with bright point sources for this rea- 
son, but even in applications with distributed emission 
regions, there can be relatively confined regions with 
larger than averaged emission (e.g. due to magnetic dis- 
sipation). In such cases, a greater degree of diffusion in 
the intensity can mitigate unphysical effects which would 
otherwise arise due to the limited angular resolution. 

A related test of an RT routine is its ability to cast a 
shadow when an optically thick obstruction is present in 
the domain. We present such a calculation in Section 5.5 
of JSD12, where the ablation of an optically thick cloud 
is studied. In this case the RT solver was used to com- 
pute the radiation field using linear interpolation for the 
intensity field of neighboring zones. Figure 14 of JSD12 
demonstrates that our RT solver can produce sharply de- 
fined umbra and penumbra under such conditions. FLD 
and other appro ximate moment method s generally fail 
this test (see e.g. lHaves fc Norma n 2003). 

5.3. Comparison with Monte Carlo and FLD Methods 

We now focus on comparing the performance of our 
short characteristics solver (referenced throughout this 
section as the SC method) with two alternative meth- 
ods: FLD and Monte Carlo (MC). Our motivation is 
two- fold: in part, we want to evaluate the performance 
on a fully three dimensional domain, but impose as few 
restrictive assumptions (e.g. the Eddington approxima- 
tion) on the radiation field as possible. Since there is a 
paucity of such truly three dimensional problems with 
analytic solutions, comparison with alternative RT solu- 
tion methods is the best alternative. In addition, FLD 
and MC methods are, in principle, some of the most com- 
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Fig. 7. — Comparison of Eddington factors computed using our RT scheme with MC and FLD computations in a representative two- 
dimensional slice near the top of a three-dimensional domain. We plot the f zz component of the Eddington tensor for computations with 
our Athena solver using either 24 (bottom, left) or 168 (bottom, center) angle bins, and for the MC (bottom, right) and FLD methods 
(top, left). We also plot p in the same slice for comparison (top, right). Note the larger range for the color bar in the panel showing the 
Eddington factor for the FLD computation. 



putationally efficient alternatives to short characteristics 
solvers, so direct comparison may allow us to assess the 
relative merits of different methods. 

For this comparison we use a three dimensional snap- 
shot from a stratified shearing box simulation, corre- 
sponding to a gas pressure d ominated patch of an accre- 
tion disk (jHirose et al.ll2006l ). This simulation was com- 
puted with t he Zeus MHD code, usi ng the FLD solver 
developed b y iTurner fc Ston e] (120 01.1 and subsequently 
modified by iHirose et al.l (|2006f) . They solved the radi- 
ation momen t equations using a flux limiter of the type 
described in iLevermore fe Pomranind (|1981[ ). Further 
details a bout the particular snapshot used here can be 
found in iBlaes et al.l ((2006). From the i? ra d dump, we 
compute -Frad and the Eddington factor f, using finite 
differences and flux limiters consistent with those em- 
ployed in the numerical simulation. 

We solve the RT equation on this snapsho t using both 
our SC solver and the MC code described in lDavis et al.l 
(2009). For both calculations, we assume isotropic elec- 
tron scattering and monochromatic RT (i.e. a single 
frequency bin) with mean opacities equal to those used 

in the Zeus simulation (x abs = 3.7 x 10 53 p n/2 E 7 e is and 
X sc = 0.33p, both in cgs units). We assume no incom- 
ing intensity at the surface boundaries and periodicity in 
the horizontal directions. The latter assumption is in- 
consistent with the use of shearing periodic boundaries 
in the radial direction in the numerical simulation, but 
this does not contribute significantly to the discrepancy 
between the SC and FLD methods 5 
We compare the radiation moments (P ra d, f, H ra d, and 

5 We have implemented shearing boundaries in our SC solver 
and confirmed this. We show the results from the SC solver with 
periodic boundary conditions to facilitate comparison with the MC 
calculation which do not support shearing periodic boundaries. 



E T ad) output by the SC/MC solvers with those deter- 
mined by the FLD method. Independent of the variable 
used for comparison, we find reasonable agreement be- 
tween the SC and MC solvers, but discrepancies with 
the FLD approximation. For brevity we will focus on a 
single scalar quantity, f zz , since it characterizes the vari- 
ation of angular distribution of the radiation field across 
methods. 

Figure [7] shows a comparison of f zz among the various 
methods for a representative two-dimensional slice near 
the top boundary of the simulation domain. In the top 
row, the left and middle panels show results from the SC 
solver, using 24 and 168 angles, respectively. The top 
right panel shows the MC results and the bottom left 
panel shows the Eddington values computed with the 
FLD approximation. The bottom right panel shows p 
for the same two-dimensional slice. 

We first compare the SC and MC calculations which 
provide similar results. The consistency of the solution 
computed by these two very different numerical meth- 
ods strongly suggests that they are providing accurate 
results. As can be seen for f zz in Figure [7] the agreement 
between the radiation moments improves as the angular 
resolution in the SC solver is increased (i.e. between the 
top left and top right panels). However, even with higher 
angular resolution there are some modest discrepancies 
in the f zz near the surface. This in part due to the statis- 
tical noise in the MC calculation, for which S/N generally 
decreases as z increases. This MC calculation was run 
with ~ 12 billion photon packets with a total compu- 
tation time that exceeded the SC solver by a factor of 
~ 100. 

Since the improvement in S/N only increases as 
roughly V~N, where N is the number of photon pack- 
ets, further improving S/N involves a substantial increase 
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in the computational time. Even for this rather large 
number of photon packets, substantial noise remains in 
the radiation field. Such a high level of statistical noise 
could lead to numerous problems when coupled to the 
MHD integrator. Hence, schemes which use MC meth- 
ods to solve RT will generally require a large number of 
packets. Our results suggest that standard MC methods 
need to be much more efficient or parallelized with ef- 
fective load balancing between the MHD integrator and 
the MC RT solver to be competitive with SC methods 
when the simulation domain is far from LTE 6 . Alterna- 
tively, it may be possible to significantly improve on this 
performance by implementing some sort of hybrid MC 
sche me to handle optically thick regions more efficiently 
(e.g. iDensmore et alj |2007) since a significant fraction of 
the time in our MC computation is spent solving RT 
in regions that are very optically thick to scattering (so 
fzz ~ 1/3) but still optically thin to absorption. 

There are several discrepancies between the SC/MC 
and FLD calculations. The most obvious is that with 
FLD, f zz approaches unity by construction in the opti- 
cally thin limit. Obtaining f zz = 1, requires the radiation 
field to be concentrated in a pencil beam of negligible 
solid angle around the z axis, and is only achieved on 
the z axis at very large distances from a finite source. 
Therefore, it is not appropriate for the upper boundary 
of a patch of an accretion disk where the radiation field is 
still rather broadly distributed over solid angle. Indeed, 
f zz ~ 0.42 is consistent with estimate s for a scatter- 
ing d ominated semi-infinite atmosphere (jChandrasekharl 
1960). In principle, one could tailor the flux-limiter to 
approach an alternative, problem dependent value, al- 
though one can imagine applications where the appro- 
priate limit will be difficult to estimate a priori. 

Furthermore, the FLD results yield f zz > 0.332 every- 
where, but in both the MC and SC calculations f zz ^ 0.3 
is frequently obtained in localized regions, consistent 
with a more horizontally directed radiation field. It is 
also clear that the FLD Eddington factors correlate with 
p to a much higher degree that in the SC or MC calcu- 
lations. Although some correlation is present in the MC 
and SC calculations as well, it is more prevalent in the 
optically thick regions and becomes much weaker in the 
optically thin regions where the radiation field should be 
more diffuse and more sensitive non-local variations in T 
and p. 

Further discrepancies between the VET and FLD ap- 
proaches are discussed in JSD12. The level at which 
these differences affect the overall dynamics and ther- 
modynamics remains unclear and ultimately requires 
comparison with full numerical simulations using the 
SC/VET methods. We note that the horizontally av- 
eraged flux in the SC and FLD methods differs by < 5% 
at the top of the domain. Hence the global thermody- 
namic properties of the simulations may not be greatly 
modified even though local properties of the radiation 
field differ. Since simulations of accretion disk dynamics 
in the shearing box approximation is one of our primary 
applications, we expect to be able t o make direct com- 
parison with FLD-based results (e.g. iHirose et alj|2006l) 

6 Although, there are problems where MC methods maybe 
preferable to SC, such as relativistic calculations that may require 
very high angular resolution if computed in the Eulerian frame. 
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Fig. 8. — Real (top) and imaginary (bottom) parts of frequen- 
cies of radiatively modified acoustic modes versus optical depth 
per wavelength to. We plot three sets of curves and symbols cor- 
responding to different Boltzmann numbers Bo = 0.01 (black), 1 
(red), 10, and 100 (blue). The curves correspond to radiation mod- 
ified acoustic modes (solid) and a non-equilibrium radiation diffu- 
sion mode (dashed). We normalize u>r and wi by the product of 
the wave number k and the adiabatic sound speed a. The real part 
of u) is zero for the radiative diffusion mode. 

in the near future. 

5.4. Radiating Linear Waves 

We now turn to tests of the RT solver when coupled 
to the MHD integrator. We first comp ute the radiative 
damp ing rate of linear (acoustic) waves ([Stein fe Spiegell 
119671 ). The closely related problem of the spatial 
damping of driven harmo nic disturbances is covered in 
iMihalas fc Mihalasl (|1984[ ). We briefly review the deriva- 
tion of the dispersion relation for such wave and refer the 
reader to these references for further discussion. We con- 
sider an ideal gas with a static, uniform background state 
in LTE, with a grey absorption opacity x and frequency 
integrated thermal source function B = ctbT^ /tt. Adopt- 
ing the notation of Mihalas fc Mihalasl ([l984f ). we define 
background and perturbed quantities with subscripts "0" 
and "1" respectively. The background states has v = 
with xo and Iq = Jq = Bq constant everywhere. 

With these assumptions and some algebra the linearly 
perturbed versions of equations CQl-fll) reduce to 



~~dt 



and 



(7-1) 



d 2 



T dpi 4vr(7 - l)xo 



Pa 9t 



Rp 



(Ji-Bi) = 0, 



dt 2 



- a?V 2 | pi 



i?p V 2 Ti = 0, 



(31) 
(32) 



where we a\ = a/y/j is the isothermal sound speed. Sim- 
ilarly, equation (jU) becomes 

n-Vh= X o{Bi-h). (33) 
To linear order we can assume 



' dB 
Bi=[ W 



Tj = 4B ^-, 



(34) 
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Fig. 9. — Real (top) and imaginary (bottom) parts of frequen- 
cies of radiatively damped linear waves versus optical d epth per 
wavelength to. The curves are analytic solutions to eq. H38H and 
the symbols are estimates derived from simulations. We plot four 
sets of curves and symbols corresponding to different Boltzmann 
numbers Bo = 0.1 (dot-dashed, squares), 1 (dashed, triangles), 10 
(dotted, crosses), and 100 (solid, diamonds). We normalize lur by 
the product of the wave number k and the adiabatic sound speed 
a. 



and solve equation (|33[) directly to evaluate J\ in equa- 
tion EH). We have 



4_Bq 
To 



7i(x 



- X0S ds, 



(35) 



where is ds is a displacement parallel to k. We consider 
plane wave solutions of the form T\ oc e J (^*- k x ) Defin- 
ing ii = k • n and integrat ing over solid angle, we obtain 
(Mih alas fc Mihalaslfl98l 



Jx 



To jo jo 
The integral evaluates to 



dfx / dy cos (kny/xo) e v . (36) 



Ji = — -tan 

T k 



(37) 



We can now solve for the dispersion relation using 
equations (JSJ), ®, and ([27]) 



uj 3 — iuj 2 vo^o ~ jafk 2 + ia 2 k 2 z/qSq = 0, 



(38) 



in ag reement with equation (16) of IStein fc Spiegell 
([19671) . We have defined 



S = 1 



■ cot 



'(f) 



and 



E. 



(39) 



(40) 



and the "0" subscript denotes that quantities are eval- 
uated using the background values. To order unity vq 
is the reciprocal of the radiative relaxation time in the 
background flow. 
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Fig. 10. — Convergence of the norm of Ll error as a function 
of resolution for linear wave in a one-dimensional domain. The 
top panel shows the error norm for waves in the adiabatic regime 
corresponding to Bo=100, to = 0.01 (crosses); Bo=100, to = 100 
(squares); and Bo=0.1, to = 10 4 (diamonds). The Bottom panels 
shows the error norm for waves in the isothermal regime corre- 
sponding to Bo=0.01, to = 0.01 (crosses) and Bo=l, to = 1 (dia- 
monds). The dashed curves show the expected trends for first-order 
(TV -1 ) and second order (N~ 2 ) convergence. 



Figure [5J shows the solutions to equation ([55)1 for var- 
ious To = Xo/k (approximately the optical depth per 
wavelength) and Boltzmann number 



Bo 



PoCpTpa _ 1670x0 
cbTo 1 v o 



(41) 



Here c p is the specific heat at constant pressure, so Bo 
is the ratio of the enthalpy flux (evaluated for v = a) to 
radiative flux. There are two types of modes: radiatively 
damped acoustic waves (solid and dashed curves) with 
phase velocity v p h = ^n/k varying between aj and a = 
y/^ai and a purely damped (ljr — 0) non-equilibrium 
radiation diffusion mode (dotted curve). 

The dimcnsionless ratio Uo'Bo/(kai) determines the im- 
portance of radiation. When this ratio is small equation 
([55)1 reduces to the standard adiabatic dispersion rela- 
tion with sound speed a and the damping rate is approx- 
imately ^0(7 — l)/(27). For ^ S /(fcai) > 1 the phase 
speed tu/k decreases, approaching the a\ when k ~ xo 
and vq 3> a\k, and the damping rate is again small com- 
pared to ka. 

For acoustic waves, the damping rate wj < ak for all Bo 
and tq. However, this is not true for radiation diffusion 



mode. For tq ~ 1, t 



i/qSq and transitions to 



t id ~ ^0=0/7 for r > 1 or t < 1. Near t ~ 1, S ~ 1, 
so the maximum decay rate has <~~> vq. If St > then 
spurious, small- amplitude oscillations may grow due to 
our failure to adequately resolve the radiative diffusion 
mode. 

Indeed, we find exactly this type of numerical instabil- 
ity for a range of To if Bo < 1. The unstable range of 
t corresponds to values for which Stc > ird f° r modes 
with wavelengths comparable to the minimum grid spac- 
ing (k ~ 1/Ax). Since we have an exact analytic solution 
for the radiation source term from equation ([37]), we can 
check this result by turning off the RT solver and up- 
dating the total energy using the exact expression for 
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Fig. 11. — Convergence of the norm of LI error as a function of 
resolution for non-grid-aligned linear wave in a three-dimensional 
domain. The crosses represent waves in the isothermal regime 
(Bo=l, to = 1) and the diamonds in the adiabatic regime (Bo=100, 
to = 0.01). As in the one-dimensional case, waves in the adia- 
batic regime converge at nearly second order, but the isothermal 
waves are nearly second order at low TV, then plateau, and finally 
transition to a regime of first-order convergence. The waves are 
computed on a 27V X N X N domain, as described in the text. 



Q ra( j. Even when the exact expression is used the code 
is numerically unstable, as expected from the argument 
above. Limiting the time step to be less than or equal to 



1 



mm — 



1 - 



cot 



(42) 

stabilizes the solution when either the exact analytic ex- 
pression or the full numerical RT solution is used to com- 
pute Q ra d. 

Since v cx x this constraint is most stringent where 
XiAxi ~ 1 in which case <5i r d = minfd/V,). Assuming 
Stc — min(Acci/ai), this implies that 



Stld ■ ra \ 
— — oc mm mo . 

otc 



(43) 



Hence, whenever the Bo number in any gridzone of the 
domain is less than unity, the maximum allowed time 
step will be determined by the radiation constraint, un- 
less some other physics (e.g. microphysical dissipation or 
magnetic fields) enforces a shorter time scale. 

We now use these solutions to evaluate the convergence 
properties of the MHD integrator when our RT solver 
is used. We simulate periodic domains with different 
combinations of tq and Bo. We initialize the background 
with v = 0, po = 1 and 7 = 5/3. The initial perturbation 
is an eigenfunction with dimensionless amplitude A — 
10~ 6 . We simulate for one adiabatic crossing time tf = 
L/a and fit for the decay rate and phase velocity. 

Figure [5] shows a comparison of the numerically de- 
rived dispersion relation with solutions of equation (|38p . 
Each symbol corresponds to fits to a simulation of a one- 
dimensional domain with TV = 256. Each curve corre- 
sponds to a different choice of Bo. We find good agree- 
ment with theory for the phase velocities OJn/k and prop- 
erly capture the transition from adiabatic to isothermal 
and back to adiabatic as to increases. The agreement 
for decay rates (~ uj 1 ) is also good except for very low 



or very high To and high Bo. In this case, the damping 
rate utj 1 is very long compared to a wave period and 
higher resolution is required to reduce the damping from 
numerical diffusion. 

We now examine convergence properties in the charac- 
teristic regimes. Figure ITUl shows the convergence of the 
norm of the LI error vector, defined as 



Sq = 



1 Ei 

N ^ 



(44) 



where q° is the eigenfunction used to initialize the do- 
main at t = to, but evaluated at t = tf. Each curve 
in Figure [TU] corresponds to a set of simulations with 
different combination of Bo and To . The plotted simula- 
tions were run on one-dimensional domains with = 4, 
but we obtain nearly identical results for grid aligned 
waves in two-dimensional (N x N) and three dimensional 
(N x N x N) domains. 

Comparison with Figure [9] shows that all of the simula- 
tions in the top panel are in the nearly adiabatic regime 
and those in the bottom panel are in the nearly isother- 
mal regime. Since radiation has only a small damping 
effect in the adiabatic regime, convergence is nearly sec- 
ond order, as when radiation is entirely absent. In the 
isothermal regime, convergence is closer to second order 
at lower resolution, but transitions to first order as reso- 
lution increases. Since we use an operator split update of 
the energy equation, first order convergence is expected 
when RT has a significant effect on the thermodynam- 
ics. Indeed, convergence is consistent with first order 
when the time step is set solely by the CFL condition 
(N > 256). For N < 128, St = 6t ld < St c and the radi- 
ation diffusion constraint sets the timestep. In this case 
St is only very weakly dependent on N. 

We also considered the convergence of non-grid- 
aligned waves in two and three dimensions. The three- 
di mensional case is n e arly i dentical to the test presented 
in lGardiner fc Stone! (|2008l) . We use a 2N x N x N pe- 
riodic domain, initialized with with a one-dimensional 
wave that has been rotated with sin a = 2/3 and sin /3 = 
2/V5 (see iGardiner fc Stond [20081 for further details). 
As in the one-dimensional case, the initial wave is an 
eigenmode with amplitude A = 10~ 6 and we use = 4. 
We again evolve the domain for one adiabatic sound 
crossing time and evaluate the Ll-error norm via 



Sq 



1 



2iV 3 



Ei 



(45) 



The convergence of the LI error as a function of N is 
shown for two waves in Figure 1111 The solid and dotted 
curves show the convergence for waves in the isothermal 
(Bo=l, t = 1) and adiabatic regimes (Bo=100, t = 
0.01), respectively. Comparison with Figure [TUl shows 
that the convergence properties are consistent with the 
one-dimensional/grid- aligned calculations. 

Further linear wave tests are presented in JSD12, al- 
though these assume the Eddington approximation and 
do not make use of the RT solver employed here. Since 
they solve the mixed frame moment equations, the char- 
acter of their numerically and analytically derived disper- 
sion relations differs from those presented here, although 
they agree qualitatively in the appropriate limit. 
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Fig. 12. — Profiles of the gas temperature (top), radiation temperature (middle), and density (bottom) versus distance for a radiation 
modified shock with Mo = 1. 2. The density, gas temp erature, and velocity are initialized with a semi-analytic planar shock solution 
computed using the methods of Lowric & Edwards (2008) (shown as continuous cu rves in each plot) while the initial radiation temperature 
is computed by the RT solver. Quantities are non-dimensionalized as described in Lowric & Edwards (2008) and discussed in section [5.51 
The red crosses indicate the initial conditions for all variables: the fact that the radiation temperature computed by the RT solver initially 
agrees with the semi-analytic solution confirms the accuracy of the RT solver. The blue crosses show the state of the variables after evolving 
the shock for t = L/ag. There is a small drift in the numerical solution due to the neglect of radiation pressure. 



5.5. Radiative Shocks 

We now consider the ability of the RT solver to model 
shocks in the presence of radiation. The physics of ra- 
diat ive shocks has been explo red by a number of authors 
(see lMihalas fc Mihalasll984l and references therein) and 
is generally well understood. However, radiating shocks 
are sufficiently complicated that simple analytic solutions 
for rad iative shocks are generall y not available. Fortu- 
nately, lLowrie fc Edwards! (|2008D (hereafter LE08) have 
developed fairly simple, semi-analytic methods for con- 
structing one dimensional planar solutions of radiating 
shocks, which are suitable for our purposes. 

LE08 construct their solutions using a grey non- 
equilibrium diffusion model of radiation hydrodynamics. 
Their treatment differs from ours in a few important 
ways. Rather than solving the RT equation © directly, 
they solve the radiation moment equations with Edding- 
ton approximation and assuming a diffusion relation for 
the radiative flux. They retain a number of velocity de- 
pendent terms which are absent in our treatment and 
include the radiation source term in the material momen- 
tum equation (our eq. This allows them to explore 
the radiation pressure dominated, which is not accessible 
with the methods discussed here (see, however, JSD12). 
Hence, our comparisons will be restricted to shock solu- 
tions with a low ratio of radiation to gas pressure and 
modest Mach numbers. 

LE08 solve a non-dimensionalized systems of equations 
with solutions that can be uniquely specified in terms of 
7, a a , Vq, k, and Mq using their notation. Here o~ a 



is the non-dimensional absorption cross section, Vq is 
roughly the ratio of radiation to gas thermal energy in the 
upstream flow, k is non-dimensional photon diffusivity, 
and Mo = v/ao is the upstream Mach number. The 
subscript "0" refers to upstream values in their notation. 

Following LE08, we examine solutions with 7 = 5/3, 
a a = 10 6 , To = 10~ 4 , and K = 1. In our notation, these 
parameters correspond to E rat \ = 10/9 x 10~ 4 -E gas , a = 

1/V3x l(T 3 c, and x tot = X abs = 1/V3 x lO^L -1 . Here, 
L is an arbitrary reference length scale and all variables 
are evaluated using their asymptotic upstream values. 

We construct one-dimensional planar shock solutions 
following the procedures outlined in LE08 and use the 
resulting profiles of p, v, and E glls to initialize our one- 
dimensional simulation domains. Since x tot L *C 1, we 
only simulate the region within a few photon mean-free- 
paths (A mfp - l/x abs ^ 0.003L) of the shock front. Since 
the semi-analytic solutions rely on the Eddington ap- 
proximation, we set = 2 (i.e. two-stream approxima- 
tion) for consistency. The radiation field at the bound- 
aries is fixed and assumes that the incoming radiation 
is in thermodynamic equilibrium with appropriate up- 
stream and downstream asymptotic temperature. We 
evolve the simulations for a time At = L/a, which is 
typically a factor of hundred (~ L/A m f p ) larger than the 
sound crossing time of the simulation domain. 

Figures [T2l and [T"3l show characteristic results for A4o = 
1.2 and 5, respectively. We use 128 and 1024 gridzoncs 
for the simulation with Aio = 1.2 and 5, respectively. We 
use a larger number for the Mq = 5 simulation to resolve 
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Fig. 13. — Same as Figure IT2l but for Mq = 5. The inset in the top panel shows a close-up view of the gas temperature near the 
shock front at t = L/ao as well as the initial semi-analytic model shifted to the right by an amount 0.0073L (the amount the shock front 
shifts in this timespan) shown as a dotted line. The shock shifts further since radiation pressure is more important at this Mach number, 
nevertheless the Zeldovich spike (and the overall shock profile) persists. 



the narrow Zeldovich spike (|Zei'Dovich fc Raizerlll967t ). 
We plot the gas temperature T gas = a 2 /{^R), radiation 
temperature T rac i = (-Erad/fli?) 1 ^ 4 , and p using the non- 
dimensional units of LE08. Here, an is the standard radi- 
ation constant and -E ra d is computed using the RT solver. 
The fact that T ra <i computed using the RT solver in the 
initial shock profile agrees well with the semi-analytic 
solutions is already an important test of our method. 
Sufficiently far upstream or downstream of the shock, 
the gas and radiation are in thermodynamic equilibrium 
with T gas = T ra d. Near the shock front, the temperatures 
deviate, with a radiation precursor upstream of the shock 
and Zeldovich spikes appearing downstream of the shock 
for the higher Mach number solutions. 

For each plot, we show two sets of curves correspond- 
ing to the initial and final profiles. Since we have initial- 
ized the simulation with stationary solutions computed 
in the shock rest frame, the material properties should 
not evolve with time. However, since our system of equa- 
tions differ from those used by LE08 to derive their so- 
lutions (in particular, we ignore the radiation pressure), 
our numerical solutions are only approximately station- 
ary. The effects of the terms we have neglected are small 
for the chosen parameters. Nevertheless, there is a slow 
but steady drift of the shock location in the downstream 
direction due to the neglect of the radiation force in the 
upstream direction. As the Mach number of the flow 
increases, the radiation force becomes increasingly im- 
portant and the shock front moves more rapidly in this 
frame. Even though the position of the shock drifts, the 
profile changes very little as the radiation source term in 
the energy equation is still well approximated. 

Further tests of radiative shocks are presented in sec- 



tions 5.2 and 5.3 of JSD12, including calculations that 
use the RT solver to compute the VET. 

5.6. Performance 

The added computational cost of using the RT solver 
is determined by a number of factors and will generally 
be problem dependent. A useful starting point is a com- 
parison of the computational cost to integrate the MHD 
equations for one timestep with the cost to perform a 
single iteration of the RT solver for a single frequency 
when run on a single processor. For a three dimensional 
domain with n M = 4 (i.e. 24 total rays) the RT solver re- 
quires ~ 40% as many operations as the CTU integrator. 
This is essentially the simplest type of problem that is of 
practical importance: an LTE grey problem with fixed 
intensity on the boundaries and an angular discretization 
that can yield a result beyond the Eddington approxima- 
tion. 

Many problems of interest will be more costly than this 
because we will need multiple iterations, multiple fre- 
quency bins, or higher angular resolution. The total cost 
of the RT solution scales approximately linearly with the 
number of frequency groups, total number of angles, or 
number of iterations, all of which are problem dependent. 
Even for LTE problems, periodic boundaries or domain 
decomposition may require multiple iterations. For most 
problems iteration will continue until the relative change 
in S (or J) is below some prescribed threshold S c and 
the total number of iterations may fluctuate from one 
timestep to the next, depending on conditions. 

For the code tests considered here, which used 5 C = 
10~ 5 , the number of iterations per timestep was < 5, 
depending on the problem, with 1-3 iterations being typ- 
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ical. Most of the tests were LTE and iteration was only 
used to handle boundary conditions. An exception is the 
uniform non-LTE atmosphere tests that were run with 
> 1000 iterations in order to obtain convergence of the 
absolute error. 

We emphasize that Figure |3] is not indicative of the 
typical number of iterations that need to be performed 
per timestep, even in highly non-LTE domains. The key 
point is that this calculation starts from an initial con- 
dition that assumes an LTE radiation field everywhere, 
even though the solution at the surface is far from LTE. 
As discussed in TF95, the main problem with the ALI 
methods used here is that they have a rather small spec- 
tral radius. Effectively, this means that it takes a rather 
large number of iterations for errors in the initial con- 
dition that span many gridzones to diminish. Since we 
are computing RT on each timestep, we already have an 
initial guess that is a reasonable approximation to the 
correct non-LTE solution. In particular, large (i.e. do- 
main scale) variations in the radiation field are usually 
already well accounted for by the solution from the pre- 
vious timestep. 

We anticipate that our initial solution of the radiation 
field before the first timestep may require hundreds to 
thousands of iterations for highly non-LTE problems (i.e. 
those with a significant fraction of zones having e <C 1), 
but that subsequent timesteps will only require a mod- 
est number (< 5) of iterations to obtain relative conver- 
gence lAS'l/S' < 10~ 3 . Our initial work on shearing box 
simulations (not reported here) supports this expecta- 
tion, although the number of iterations depends some- 
what on just how non-LTE the radiation field becomes. 
lHayek et al.l (|2010t ) report similar numbers of iterations 
(see their Figure 2) as being typical of their scattering 
dominated calculations. 

A second consideration affecting performance is the 
maximum timestep that can be used with the operator 
split update of the total energy described in section 01 
The generalized CFL conditions derived in section 15.41 
may reduce the timestep when E ra< \ is a significant frac- 
tion of -Egas or v <C c. For such problems it will be more 
efficient to use the VET method of JSD12 when feasible. 

A third consideration is scalability. Since scaling effi- 
ciency will be somewhat machine dependent, we are pri- 
marily concerned with assessing how the code performs 
with RT relative to the (M)HD configuration with no so- 
lution of RT. To make the comparison concrete, we study 
the weak scaling for 32 3 gridzones per core on the SciNet 
General Purpose Cluster 7 , which consists of eight core 
nodes made from two 2.53 GHz quad-core Intel Xeon 
5500 Nehalem processors. Tests were performed using 
the Infiniband interconnect. We use a grid aligned ra- 
diating linear wave on a three-dimensional domain with 
rip = 4, To = 1, and Bo=l in the radiating case and a 
slow magnetosonic wave for MHD only calculation. Both 
tests were performed using the full MHD CTU integra- 
tor. We initialize the radiating wave as described in sec- 
tion 15.41 but we only allow a single iteration of the RT 
solver per timestep. We compute the efficiency by di- 
viding the number of zone cycles/second obtained for a 
problem run with multiple cores by the number of zone 
cycles/second for a single core. The resulting scaling with 

7 http://www.scinet.utoronto.ca 



number of cores is nearly identical in the radiating and 
non-radiating cases, falling to about 75% 8 at 512 cores 
(64 nodes). 

We also find the same scaling efficiency for n M = 12, 
which is notable because this correspond to a factor of 
seven increase in the number of specific intensity bins 
that need to be passed on each iteration. In this case, 
the increase in communication demands is balanced by 
the roughly factor of seven increase in the cost of comput- 
ing the RT solution with more angle bins. These results 
suggest that the scaling efficiency with the RT solver 
will tend to follow the non-radiative scaling when they 
are run on the same platform, as long as the number of 
iterations remains constant. 

The assumption that number of iterations stays fixed 
is an important caveat. In fact, this assumption will 
not generally hold since we use iteration to handle sub- 
domain boundaries. To understand this, it is useful to 
consider two LTE problems: one in which the optical 
depth across each subdomain is very large and one for 
which the whole domain is optically thin. In first prob- 
lem, the radiation field is determined entirely locally and 
propagation of changes in the radiation field from neigh- 
boring subdomains will require only a single iteration. In 
the second problem, variations in the emissivity on one 
side of the domain can modify the radiation field on the 
other. For a cubic array of N subdomains, it could take 
~ TV 1 / 3 iterations to propagate the radiation across the 
entire domain. This means that scaling efficiency can, in 
principle, be very problem dependent. For most of the 
applications of primary interest, the majority of subdo- 
mains will be optically thick so increasing the number 
of subdomains should not significantly increase the num- 
ber of iterations required. Therefore, scaling efficiency 
should reasonably consistent with the non-radiating case. 

Performance and scaling of the overall VET scheme is 
discussed in section 7.2 of JSD12. 

6. SUMMARY 

We have described our implementation of an RT solver 
in the Athena MHD code. Our module implements a 
short characteristic method for computing RT on Carte- 
sian, multidimensional simulation domains. The RT 
equation is solved once each simulation timestep for a 
computational cost comparable to or less than a single 
timestep of the MHD integrator for simple (e.g. LTE 
grey) problems. Since we are focused on astrophysical 
problems where velocities are slow compared with the 
speed of light, we drop the time derivative of intensity 
and the system becomes an integrodiffcrential equation 
with no explicit time dependence for the radiation field. 
The material properties of the flow are evolved using the 
standard Athena MHD integrators, but with radiative 
heating and cooling source terms computed from the RT 
solver. The code solves the RT equation for frequency 
dependent, absorption and scattering opacities. Non- 
LTE effects arising from scattering processes are han- 
dled with ALI methods. The resulting code is well-suited 
for non-relativistic astrophysical problems where diffuse 
emission, rather than bright point sources, constitutes 

8 Note that this scaling efficiency for the non-radiating wave is 
somewhat lower than previous tests on other platforms. See e.g. 
https: / /trac. princeton.edu/Athena/wiki/AthenaDocsScaling. 
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the dominant source of radiation. 

We provide a detailed summary of the short charac- 
teristics and ALI implementations. We also describe our 
implementation of an operator split update of the en- 
ergy equation using a radiation source term computed 
directly by the RT solver. Alternatively, the RT solver 
can be used to compute a VET, which can then be input 
into the integration of the coupled MHD and radiation 
moment equations. The use of the RT module for this 
purpose is discussed in a companion paper (JSD12). 

We also present results from several test problems, 
which roughly fall into two classes: tests of the RT solver 
on static simulation domains, and tests of the coupled 
RT solver and MHD integrator for time-dependent hy- 
drodynamics simulations. These tests demonstrate the 
accuracy of the RT solver for multidimensional problems, 
assess its convergence properties for applications where 
scattering leads to significant deviations from LTE. They 
indicate that substantial improvements in accuracy and 
efficiency may be obtained over alternative methods, 
such as flux-limited diffusion and Monte Carlo based RT 
solvers. 

The tests also evaluate the accuracy and stability of 
the MHD integrator when coupled to the RT solver via 
operator splitting. They verify that the code is generally 
only first-order accurate for problems where heating or 
cooling of the fluid by the radiation field is significant. 
They also illuminate important time step constraints and 
are useful for assessing the efficiency and accuracy of the 
code for solving various astrophysical problems. 

In particular, we derive a generalized CFL condition, 
predicated on the need to resolve the non-equilibrium ra- 
diation diffusion mode. The requirement to place some 
limi ts on timeste ps due to rapid radiative relaxation (see 
e.g. [Castor 2003) are generally acknowledged and imple- 
mented in previous work. However, we have not seen any 
explicit reference to resolving the damping rate of the ra- 
diation diffusion mode, which provides a a practical and 
precise guideline for ensuring numerical stability. 

The focus in this work has been on modeling the fluid 
dynamics and thermodynamics with a self-consistent 
computation of the radiation field. However, we ex- 
pect that the ability of the RT solver to provide de- 
tailed outputs of the emergent radiation field, such as 
images, lightcurves, and spectra, may be equally impor- 
tant. Indeed, for some applications the production of 
such outputs may be the primary motivation for includ- 

9 http:/ /trac. princeton.edu/Athena 



ing radiation in the simulation. In principle, the RT 
solver can be used solely to generate diagnostic outputs, 
even in simulations without self-consistent feedback of 
radiation on the material flow, either in real time or via 
post- processing. 

In addition to the simple test problems described here, 
we are beginning a research program to simulate the 
local structure of accretion flows (i.e. stratified shear- 
ing boxes) with radiative heating and cooling. Appli- 
cations to radiation dominated environments using the 
VET method (JSD12) are also underway, and include 
studies of radiation dominated accretion disks, radiative 
Ray leigh- Taylor instability, and the radiative driving of 
cold gas. Applications to extrasolar planets, star form- 
ing environments, boundary layers, galactic and accre- 
tion disk outflow are also under consideration for future 
work. 

A significant limitation of the RT solver described here 
is the short-characteristics method's inability to accu- 
rately handle bright points sources. We plan to address 
this in future work using a hybrid scheme that computes 
the direct radiatio n from point sou r ces usi ng the algo- 
rithms described in lKrumholz et all ()2007bl ) and models 
the diffuse emission with RT solver described here. 

The source code for our RT solver, test problems and 
associated documentation will be included in the publicly 
available version of Athena 9 in the near future. We have 
endeavored to make the code as user friendly as possible 
and strongly encourage interested parties to use the code 
in their own research. 
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